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Abstract 

We  consider  the  qualitative  behaviour  of  exact  and  approximate  solutions  of  integral 
and  integro-differential  equations  with  fading  memory  kernels.  Over  long  time  intervals 
the  errors  in  numerical  schemes  may  become  so  large  that  they  mask  some  important 
properties  of  the  solution.  One  frequently  appeals  to  stability  theory  to  address  this 
weakness,  but  it  turns  out  that,  in  some  of  the  model  equations  we  have  considered, 
there  remains  a gap  in  the  analysis. 

We  consider  a linear  problem  of  the  form 

y'(t)  - - f e~x{t~s)y(s)ds,  2/(0 ) = 1, 
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and  we  solve  the  equation  using  simple  numerical  schemes.  We  outline  the  known  sta- 
bility behaviour  of  the  problem  and  derive  the  values  of  A at  which  the  true  solution 
bifurcates.  We  give  the  corresponding  analysis  for  the  discrete  schemes  and  highlight 
that,  for  particular  stepsizes,  the  methods  give  unexpected  behaviour  and  we  show  that, 
as  the  step  size  of  the  numerical  scheme  decreases,  the  bifurcation  points  tend  towards 
those  of  the  continuous  scheme.  We  illustrate  our  results  with  some  numerical  examples. 


1 Introduction 

The  qualitative  behaviour  of  numerical  approximations  to  solutions  of  functional  differ- 
ential equations  is  an  important  area  for  analysis.  We  aim  to  investigate  whether  the 
behaviour  of  the  numerical  solution  reflects  accurately  that  of  the  true  solution.  We  are 
particularly  concerned  with  the  behaviour  of  the  solution  over  long  time  periods  when  (in 
particular)  the  convergence  order  of  the  method  gives  us  limited  insight,  since  the  error 
depends  on  a constant  that  grows  with  the  time  interval.  Many  authors  are  concerned 
with  stability  of  solutions  and  of  their  numerical  approximations.  We  have  considered 
elsewhere  (see  [7])  the  stability  of  numerical  solutions  of  equations  of  this  type  (and  of 
non-linear  extensions).  This  analysis  raised  a number  of  questions,  which  we  consider 
here,  about  just  how  well  the  full  range  of  qualitative  behaviour  of  even  quite  a simple 
equation  is  understood. 

Bifurcations  (by  which  we  shall  mean  any  change  in  the  qualitative  behaviour  of 
solutions)  frequently  arise  only  for  systems  or  for  higher  order  problems  and  therefore 
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one  is  particularly  interested  in  finding  suitable  simple  equations  as  the  basis  for  analysis. 
In  this  paper,  we  consider  the  solution  by  numerical  techniques  of  the  integro-differential 
equation 

•1) 

The  equation  is  a linear  convolution  equation  with  separable  fading  memory  con- 
volution kernel  and  therefore  is  a simple  example  from  an  important  class  of  problems 
familiar  in  applications.  It  is  also  possible  to  analyse  the  equation  in  the  form  of  a second 
order  ordinary  differential  equation. 

The  equation  has  several  key  properties  that  make  it  an  ideal  basis  for  our  analysis: 

1.  it  depends  on  the  value  of  the  single  parameter  A, 

2.  when  A varies  through  real  values,  four  distinctive  qualitative  behaviours  in  the 

solution  can  be  detected,  and 

3.  equations  with  exponential  convolution  kernels  frequently  arise  in  applications  and 

elsewhere  in  the  literature. 

For  A real  and  positive,  the  kernel  is  of  fading  memory  type.  For  A real  and  negative, 
the  kernel  has  a growing  memory  effect.  This  linear  equation  displays  surprisingly  rich 
dynamical  behaviour  for  real  values  of  the  parameter  A and  it  is  this  behaviour  that  we 
want  to  consider  for  the  numerical  scheme.  We  note  that  the  classical  test  equation 

y'(t)  = g{t)  + &/(f)  + 1?  f y(s)ds , rj  ^ 0 (1.2) 
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([1,  2])  displays  the  same  range  of  qualitative  behaviour  possibilities  as  (1.1)  for  varying 
values  of  the  two  real  parameters  £,  r/. 

This  motivates  us  to  consider  equation  (1.1)  as  a prototype  problem  that  is  interesting 
in  its  own  right  and  that  will  also  provide  insight  into  the  behaviour  of  more  complicated 
equations.  We  propose  to  give  a further  analysis,  where  we  consider  the  boundaries  along 
which  bifurcations  occur  for  equation  (1.2)  in  a sequel  [3]. 

We  consider  the  following  questions. 

1.  Does  the  numerical  scheme  display  the  same  four  qualitatively  different  types  of 
long  term  behaviour  as  are  found  in  the  true  solution? 

2.  Are  the  interval  ranges  for  the  parameter  A that  give  rise  to  the  changes  in  behaviour 
of  the  solution  the  same  as  in  the  original  problem? 

2 Behaviour  of  the  exact  solution 

We  consider  the  equation  (1.1)  which  can  be  shown  to  have  a unique  continuous  solu- 
tion (see,  for  example,  [10]).  One  can  easily  establish  (by  considering,  for  example,  an 
equivalent  ordinary  differential  equation)  the  general  solution 

. . -A+\/a2-1,  -A-v'A2-4j  , 

y(t)  = Ae * t + Be 2 * (2.1) 

where  A , B are  constants.  For  real  values  of  A the  solution  to  (1.1)  bifurcates  (or  changes 
qualitative  behaviour)  at  A = 0,  ±2.  We  have  the  following  qualitative  behaviour. 


/(*)  = - / 
Jo 


e A(t  S)y{s)ds,  y{ 0)  = 1. 
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Al,  When  A > 2,  y —>  0 as  t — too,  with  no  oscillations. 

A2.  When  0<A<2,  y — ► 0 as  t — > oo,  with  infinitely  many  oscillations. 

A3.  When  A = 0,  y(t)  = cos (t)  (persistent  oscillations). 

A4.  When  — 2 < A < 0,  the  solutions  contain  infinitely  many  oscillations  of  increasing 
amplitude. 

A5.  When  A < —2,  the  solution  grows  (in  magnitude)  without  any  oscillations. 

3 Numerical  analysis 

To  apply  a numerical  method  to  an  integro-differential  equation  of  the  type 

y'(t)  = f (t,y(t),  k{t,s,y{s))dsj  , y(0)  = yo,  (3.1) 

we  write  the  problem  in  the  form 

y'(t)  = f{t,y(t),z(t))  (3.2) 

z{t)  — I k(t,s,y(s))ds.  (3.3) 
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We  solve  (3.2),  (3.3)  numerically  using  a linear  multistep  method  for  solving  equation 
(3.2)  combined  with  a suitable  quadrature  rule  for  deriving  approximate  values  of  z from 
equation  (3.3)  (see  [2]).  Such  a method  is  sometimes  known  as  a DQ-method.  For  linear 
fc-step  methods,  one  also  needs  to  provide  a special  starting  procedure  to  generate  the 
additional  k — 1 initial  approximations  to  the  solution  that  are  not  given  in  the  equation 
but  are  needed  by  the  multistep  method  on  its  first  application.  It  turns  out  that  one 
needs  to  choose  the  quadrature,  multi-step  method  and  starting  schemes  carefully  to 
ensure  that  the  resulting  method  is  of  an  appropriate  order  of  accuracy  for  the  work 
involved.  One  should  try  to  choose  schemes  of  the  same  orders  as  one  another  since  the 
order  of  the  overall  method  is  equal  to  the  lowest  of  the  orders  of  the  three  separate 
methods  (the  multistep  formula,  the  starting  value  scheme  and  the  quadrature)  used 
to  construct  it.  In  this  paper  we  have  chosen  to  focus  on  one-step  methods.  There  are 
two  reasons  for  this:  we  have  thereby  avoided  the  need  to  construct  special  starting 
procedures  which  would  make  our  analysis  more  complicated;  as  Wolkenfelt  showed  in 
[11],  methods  with  a repetition  factor  of  1 (such  as  the  ones  we  consider)  are  always 
stable  and  we  also  draw  attention  (see  [9]  for  example),  to  the  fact  that  the  trapezoidal 
rule  is  an  A-stable  1-step  method. 

For  a well-behaved  numerical  scheme  for  (3.2),  (3.3),  we  would  anticipate  four  in- 
tervals (as  with  the  continuous  problem)  of  A-values  where  the  solutions  to  the  discrete 
scheme  behave  qualitatively  differently.  However  we  know  from  investigation  of  bifurc- 
ation points  for  numerical  solution  of  delay  differential  equations  (see  [12])  and  indeed 
from  stability  analysis  of  integro-differential  equations,  that  the  points  at  which  the  qual- 
itative behaviour  of  the  solution  changes  may  arise  at  the  wrong  values  of  the  parameter. 
Based  on  previous  experience  (see  [6])  we  would  expect  this  difference  to  be  dependent 
upon  the  stepsize  h of  the  numerical  method  and  on  the  choice  of  method  itself.  Further- 
more (see,  for  example  [8],  [12]),  one  might  expect  the  bifurcation  points  of  the  discrete 
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scheme  to  approach  the  bifurcation  points  of  the  continuous  problem  as  h — > 0 and  one 
could  anticipate  that,  for  a method  of  overall  order  p,  the  approximation  of  the  true  bi- 
furcation point  by  the  bifurcation  point  of  the  numerical  scheme  would  also  be  to  0(hp). 
We  will  show  in  this  paper  that  (for  h — » 0)  the  approximation  of  the  bifurcation  points 
in  the  methods  we  have  chosen  is  at  least  to  the  order  of  the  method. 

To  keep  the  analysis  reasonably  simple,  we  consider  the  following  discrete  form  of 
(3.2).  We  use  a linear  0-method  in  each  case  so  that  we  solve  the  system 

Vn+l 

Fn 

Zn 

One  could  choose  any  combination  of  0j,O  < 0i  < 1 and  a natural  choice  could  be 
0i  =02-  However,  in  order  to  start  with  a simple  method  where  the  algebraic  problem 
is  tractable  we  have  considered  first  the  cases  where  0i  = 0 and  we  consider  a range  of 
values  of  02. 

One  solves  equations  of  the  form 

2/71+1  - Vn  = ~h2  ^e-Hn+Vhy0  + ^e~^(n+l 

Note  that  we  have  used  a simple  procedure  to  find  the  additional  starting  value  yi  = 1. 
We  have  observed  from  the  integro-differential  equation  that  y'{ 0)  = 0 and  have  deduced 
that  y(h ) = y( 0)  will  provide  a reasonable  order  1 starting  approximation.  This  choice 
of  formula  implies  that  we  are  combining  a backward  Euler  scheme  to  discretise  the 
differential  equation,  with,  respectively,  (for  02  = 1)  the  forward  rectangular  (Euler) 
rule,  (for  02  = |)  the  trapezoidal  rule  and  (for  02  = 0)  the  backward  rectangular  rule 
for  the  quadrature.  We  will  return  to  consider  other  combinations  of  0i,02  later. 

The  equation  (3.7)  is  equivalent  to 

(1  + h2(  1 - 02))  t/n+2  + ( h262er>'h  - 1 - e~Xh)  yn+1  + e~Xhyn  = 0.  (3.8) 

The  behaviour  of  the  solution  as  t — > oo  depends  on  the  roots  of  the  characteristic 
equation 

(l  + h2(l  - 02))  k2  + (h202e-Ah  — 1 — e~Xh)  k + e~Xh  = 0.  (3.9) 

Any  solution  of  (3.8)  will  be  asymptotically  stable  if  both  roots  of  (3.9)  are  of  magnitude 
less  than  one  and  unstable  if  either  root  of  (3.9)  has  magnitude  greater  than  one.  The 
solutions  will  contain  (stable  or  unstable)  oscillations  when  the  roots  of|(3.9)  are  complex 
or,  indeed,  when  at  least  one  root  is  negative.  It  follows  from  this  ■’’(see  [4])  that  the 
bifurcations  occur  as  follows  (for  reasonably  small  h > 0). 


j)yj  + (l  - 02)y„+i^  , yq  = yi  - !• 

(3-7) 


— Vn  + h(6iFn  + (1  — 9i)Fn+i),  n — 0, 1, ... , (3.4) 

= f(nh,yn,zn),  (3.5) 
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di  ,,r,  , . i,  / i+2h’2-h2e2-2J-h2(h2e2-i-h2)\  . ... 

Bl.  When  A > -ft  In  I h^e2'-2h2e2+ 1 J , y„  -+  0 as  n —>  oo  with  no 

oscillations.  This  condition  can  be  written  in  the  simpler  form 

A > i In  (l  + 2h2  - h262  + 2a/— ft2  ( h?e2  - 1 - fi2)) 

and  we  thank  the  anonymous  referee  for  pointing  out  this  simplification. 

B2.  When  ft  In  < A < ft  In  (l  + 2 h2  - h262  + 2 y/-h2  (/i202  - 1 - /i2)) , 

yn  — > 0 as  n — > oo,  with  infinitely  many  oscillations. 

B3.  When  A = ft  In  ^ i+p~(i-6>2) ) we  obtain  persistent  oscillations. 

B4.  When  ft  In  (l  + 2 h2  - h292  - 2/-/i2  (h202  - 1 - ft2))  < A < ft  In  (1+„^_e2)), 
the  solutions  contain  infinitely  many  oscillations  of  increasing  amplitude. 

B5.  When  A < ft  In  ^1  + 2 h2  — h292  - 2/— ft2  ( h?02  — 1 — ft 2)),  the  solution  grows  (in 
magnitude)  without  any  oscillations. 

4 Bifurcation  points  of  the  numerical  scheme  as  approximations 
to  true  bifurcation  points 

We  consider  now  the  way  in  which  the  bifurcation  points  of  the  discrete  scheme  approx- 
imate those  of  the  original  problem.  We  are  using  a numerical  scheme  of  order  1. 

First  we  consider  the  value  of  Ai  = ft  In  ^1  + 2 h2  - h202  + 2/ -h2  (h202  - 1 - h?)j 
as  d2  varies  and  h — > 0.  It  is  easy  to  see  that,  as  fi  — ► 0,  the  value  X\  satisfies  A]  — * 2. 
In  fact  we  can  give  greater  precision  to  this.  We  can  show  that  Ai  = 2 - 02h  + 0({€) 
as  h —*  0.  This  means  that,  for  9 methods  in  general,  the  approximation  by  our  scheme 
approximates  the  true  value  (—2)  to  order  1 (the  order  of  the  method),  as  h — > 0.  In  the 
particular  case  02  — 0 the  approximation  is  to  order  2. 

For  A2  = ft  In  ^1+//x_g2^  it  is  straightforward  to  show  that  stability  is  lost  at  a 
value  of  A that  approximates  the  true  value  (0)  to  order  1 in  general.  In  fact,  for  02  = 1, 
the  forward  Euler  scheme,  the  approximation  is  exact  for  all  values  of  h. 

The  analysis  of  A3  = ft  In  ^1  + 2 h2  — h262  — 2 /— fi2  (h292  — 1 - h2)J  follows  in  ex- 
actly the  same  way  as  for  Ai  and  leads  to  an  identical  conclusion:  the  approximation  of 
the  bifurcation  point  A = — 2 is  in  general  to  order  1 as  h — > 0 and  to  order  2 if  62  = 0. 

We  illustrate  our  results  graphically.  Each  of  the  plots  shown  in  Figure  1 illustrate, 
for  varying  h,  the  ranges  for  the  parameter  A where 

1.  the  solutions  are  unstable  due  to  at  least  one  real  root  greater  than  unity  in  magnitude 

(the  darkest  region  in  the  figures)  (exponential  growth  if  the  root  is  positive,  growing 
oscillations  if  the  root  is  negative), 

2.  the  solutions  are  unstable  due  to  growing  oscillations  (the  next  darkest  region  in  the 

figures), 

3.  the  solutions  are  stable  with  asymptotically  stable  oscillations  (the  lightest  region  in 

the  figure),  and 
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4.  the  solutions  are  stable  with  exponentially  stable  decay. 

We  can  compare  with  the  right  hand  plot  in  Figure  2 which  shows  the  true  regions 
for  the  original  problem  and  we  can  make  the  following  observations. 

1.  As  ft  — * 0 the  values  of  A at  which  changes  in  the  behaviour  occur  approach  the  true 

values.  This  coincides  with  our  previous  experience  in  delay  equations  (see  [8]). 

2.  There  is  some  extremely  surprising  behaviour  for  some  values  of  ft  > 0. 

(a)  For  the  two  values  62  = 0.5  and  62  — 1 we  can  see  that  the  darkest  region  is  in 
two  parts:  in  the  upper  part  there  is  a negative  real  root  of  magnitude  greater 
than  unity  leading  to  exponentially  growing  oscillations  in  the  solution;  in  the 
lower  part  there  is  a positive  real  root  of  modulus  greater  than  unity  leading 
to  exponential  growth  in  the  solutions. 

(b)  There  can  be  a critical  value  of  ft  > 0 (ft  = ^==  when  62  > 0)  at  which,  for 
apparently  arbitrarily  large  A < 0 the  numerical  solution  displays  oscillatory 
behaviour. 

(c)  There  can  be  an  additional  thin  region  (visible  only  in  larger  scale  versions 
of  the  plots)  between  the  darkest  and  lightest  regions  in  which  there  is  a real 
negative  root  of  magnitude  less  than  unity  leading  to  decaying  oscillations. 

(d)  For  62  = 0.5  and  62  = 1 the  upper  part  of  the  darkest  region  indicates  some 
really  strange  behaviour:  spurious  oscillations  may  arise  for  arbitrarily  large 
negative  values  of  A and  even  (see  figure  1)  for  some  positive  values  of  A.  Thus 
we  can  have  the  situation  (for  example  for  A small  and  positive)  where  the  true 
solution  tends  to  zero  while  the  approximate  solution  exhibits  oscillations  of 
growing  magnitude.  Alternatively,  (for  A large  and  negative)  the  true  solution 
could  exhibit  high  index  exponential  growth  while  the  approximate  solution 
exhibits  oscillations.  We  draw  attention  also  to  the  fact  that,  for  O2  — 0.5 
and  62  = 1 the  stability  boundary  of  the  method  is  made  up  of  parts  of  the 
boundaries  of  two  regions,  making  the  prediction  of  behaviour  for  varying  ft  > 0 
particularly  difficult. 

We  believe  that  these  observations  justify  our  view  that  more  attention  needs  to 
be  paid  to  changes  in  qualitative  behaviour  other  than  stability  in  reaching  a good 
understanding  of  the  behaviour  of  numerical  methods  for  problems  of  this  type. 

We  can  consider  next  whether  these  observations  are  equally  true  for  other  choices  of 
numerical  method.  We  present  in  Figures  2 plots  revealing  the  qualitative  behaviour  of 
solutions  to  equations  (3.2),  (3.3)  with  other  choices  of  0-method.  It  is  easy  to  see  that, 
even  for  combinations  such  as  using  the  trapezium  rule  for  both  parts  of  the  discretisation 
(a  method  characterised  by  9\  = 62  = 0.5  and  known  to  do  very  well  at  preserving  the 
stability  boundary)  there  are  problems  in  the  preservation  of  other  types  of  qualitative 
behaviour  when  ft  is  not  very  small.  Similarly,  we  can  see  that  the  choice  9\  — 62  = 1 
leads  to  a shrinking  range  (as  ft  increases)  for  A that  lead  to  stable  oscillatory  solutions. 
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Fig.  2.  Bifurcation  points  as  h varies  for.  respectively, 
analytical  problem. 
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5 Alternative  approaches 

The  particular  equation  we  have  considered  can  be  formulated  as  an  integro-differential 
equation,  as  an  integral  equation  or  as  a second  order  differential  equation.  We  have 
shown  in  [4]  that  the  interesting  and  somewhat  surprising  observations  about  numerical 
behaviour  that  we  made  in  the  previous  section  also  apply  in  these  other  formulations. 

6 Closing  remarks 

The  results  presented  in  this  paper  show  that  the  well-established  stability  theory  based 
on  the  analysis  of  equation  (1.2)  gives  only  a very  limited  insight  into  the  qualitative 
behaviour  of  solutions  of  the  class  of  convolution  equations  with  exponential  memory 
kernel  that  we  have  considered  here.  We  have  observed  elsewhere  (see  [5,  6,  7])  that  the 
qualitative  behaviour  of  numerical  solutions  to  equations  of  this  type  may  have  surprising 
features  and  our  consideration  here  of  the  prototype  problem  (1.1)  illustrates  how  this 
unexpected  behaviour  may  arise.  We  have  seen  in  this  paper  how  oscillations  may  arise 
in  the  numerical  schemes  when  they  should  not,  and  how  in  other  cases  the  numerical 
schemes  may  supress  genuine  oscillatory  behaviour.  When  one  seeks  good  methods  based 
on  a stability  analysis,  the  desire  is  to  focus  on  those  methods  where  the  step-length  h > 0 
is  not  subject  to  some  upper  bound  to  ensure  the  stability  of  the  method.  However  our 
initial  observations  in  this  paper  have  shown  that  this  may  well  prove  an  unreasonable 
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expectation  when  one  is  investigating  these  other  changes  in  qualitative  behaviour. 

We  believe  that  this  paper  introduces  a range  of  worthwhile  investigations  in  a field 
that  is  still  quite  open.  Space  restrictions  have  prevented  us  from  considering  the  be- 
haviour of  more  general  methods  in  this  paper  and  also  from  extending  our  analysis  to 
consider  other  problems.  The  results  we  have  presented  here  show  that,  for  these  simple 
methods  at  least,  the  bifurcation  parameters  are  approximated  in  the  numerical  scheme 
to  at  least  the  order  of  the  method,  for  sufficiehtly  small  h > 0.  It  is  also  very  clear 
that,  even  for  what  appears  to  be  a simple  problem,  the  choice  of  numerical  scheme  and 
the  form  in  which  the  problem  is  presented  provide  us  with  a rich  source  of  example 
behaviour. 
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